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Ordering of the frustrated classical Heisenberg model on the triangular-lattice with an incom- 
mensurate spiral structure is studied under magnetic fields by means of a mean-field analysis and 
a Monte Carlo simulation. Several types of multiple-g' states including the skyrmion-lattice state 
is observed in addition to the standard single-^ state. In contrast to the Dzyaloshinskii-Moriya 
interaction driven system, the present model allows both skyrmions and anti-skyrmions, together 
with a new thermodynamic phase where skyrmion and anti-skyrmion lattices form a domain state. 

PACS numbers: 75.10.Hk, 05.50.+q, 75.40.Mg, 64.60. F- 



Ordering of geometrically frustrated magnets has at- 
tracted recent interest [l|-[3| . The antiferromagnetic ( AF) 
Heisenberg model on the triangular lattice is a typi- 
cal example of two-dimensional (2D) geometrically frus- 
trated magnets. When AF nearest-neighbor (NN) inter- 
action only, the ground state of this system is the three- 
sublattice 120° structure, which is commensurate to the 
underlying lattice. Under magnetic fields, the ordered 
state still keeps the three-sublattice structure and leads 
to a rich phase diagram [4, 5]. 

When further- neighbor interactions become dominant, 
the ground state often takes an incommensurate spi- 
ral structure. In recent experiments on the triangular- 
lattice AF NiGa2S4, an incommensurate spiral state due 
to strong further-neighbor interaction has been reported 
jsl-lsl. In this compound, magnitude of the AF third- 
neighbor interaction | J3 1 is larger than the ferromagnetic 
NN interaction Ji, \Jz\/Ji ~ 5. In another compound 
NiBr2, an incommensurate state due to the AF third- 
neighbor and the ferromagnetic NN interactions were re- 
ported where | J3I/ Ji - 0.262 d, [io|. 

From a symmetry viewpoint, an important difference 
of the incommensurate spiral from the 120° structure is 
that the ground state possesses a three-fold degeneracy 
with respect to the choice of three equivalent directions 
of wavevectors on the lattice. Generally, a phase transi- 
tion related to the breaking of such a discrete degeneracy 
could occur even in a 2D Heisenberg model. Indeed, the 
classical J1-J3 model in zero field exhibits a first-order 
transition associated with a breaking of such three-fold 
C3 lattice symmetry. The ordered state is a single-g state 
where one of three equivalent wavevector directions is 
chosen lll-[l3|. 

This three-fold degeneracy could also be a source of 
exotic ordered states, e.^., various types of multiple-(7 
states where more than one wavevectors coexist. Al- 
though multiple-g states were not reported in previous 
zero-field calculations, it might be realized under applied 
fields. Multiple-g states are generally incompatible with 
the fixed spin- length condition \Si\ = 1 and not favored 



at lower temperatures in the classical system, whereas 
they might be stabilized at moderate temperatures by 
thermal fluctuations. 

In this letter, based on a mean-field calculation and 
a Monte Carlo simulation, we show that several types 
of multiple-g' states are indeed stabilized under magnetic 
fields, and that one of them corresponds to the so-called 
"skyrmion-lattice" state l3-23|. 

We focus here on the triangular-lattice Ji- J3 ior J1-J2 
model jin a magnetic field of intensity H whose Hamil- 
tonian is given by 



Sj-HY,Si,., (1) 



where Xl(^j) and Xl^^^^-^^ mean the sum over the NN 
and the third-neighbor (or the second-neighbor) pairs, 
respectively. We consider classical Heisenberg spins 
Si = {Si^x-, Si^y^ Si^z) with \Si\ = 1. An incommensu- 
rate ground state appears, in the J1-J3 model, for the 
ferromagnetic NN interaction Ji > and the AF third- 
neighbor interaction J3 < with J1/IJ3I < 4, while, in 
the J1-J2 model, for the AF J2 with -1 < Ji/| J2I < 3. 

We first perform a mean-field analysis based on the 
method of Reimers et al. 24, 25[. Up to quartic order. 



the Landau free energy of the model is given by 
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where is the order parameter corresponding to the 
Fourier magnetization given by = (Sq) with Sq = 
W '^i exp(— • Ti) {N the number of spins). The sum 
^1^1 runs over all gr^s satisfying ^^Qi = 0, and Jq is 
the Fourier transform of the exchange interaction. 

From the quadratic term of the free-energy expan- 
sion, one sees that the Fourier mode corresponding to 
the maximum Jq becomes unstable at Tc = \Jq* where 
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gf* is the critical wavevector. In the J1-J3 model, gr* 
appears along the direction of NN bonds with \q*\ = 

I ^1 + ^1 — ^j^^ {a is the lattice constant), 



■ cos 



while in the the J1-J2 model it appears along the second- 
neighbor direction with |gr*| = cos~^ + tJ") • 
Note that, reflecting the C3 symmetry of the lattice, gr* 
is degenerate as i^i, ig^l and ig'J. 

Just below the transition temperature, one might ne- 
glect all wavevectors other than the three critical modes 
gr* and the uniform gr = mode. Within this approxima- 
tion, we find three metastable ordered states in addition 
to the paramagnetic state. These ordered states are char- 
acterized by the number of wavevectors appearing in the 
xy component; they are single-g, double-g, and tiiple-q 
states. 

(i) The single-q state: Spins form an umbrella structure 
where the xy component forms a single-g spiral charac- 
terized by one of the three gr*s, while the z component 
consists of a uniform q = component along H. This 
state is compatible with the condition = 1. In fact, 
the ground state of our model turns out to be this single-g 
state irrespective of the field intensity. 

(ii) The douhle-q state: Here the xy component is a su- 
perposition of two spirals with, e.g, igr^ and ±gr|, while 
the z component forms a lineally polarized spin density 
wave with ±gr| complementary to the xy component with 
an additional uniform component. 

(iii) The triple-q state: This state is a superposition of 
three (distorted) spirals characterized by wavevectors ql , 
gr|, and gr| with an additional uniform component along 
H. In contrast to the single-g' and the double-g states, 
three spiral planes are perpendicular to the xy plane and 
rotate by 120 degrees with each other. By using arbitrary 
three unit vectors lying on the xy plane, {i = 1,2,3) 
satisfying = 0, it is given by 



Si^xy = Ixy ^ sin(gr* • n + Oj)ej, (3a) 
3 

Si,z = cos(gr* • + 6>^) + m^, (3b) 



where we introduced the phase factors 6i {i = 1 ^ 3) 
satisfying the constraint cos(^i + 6^2 + ^3) = — 1, Ixy and 
Iz being T-dependent constants. 

Interestingly, this triple-g state spin configuration just 
corresponds to the ''skyrmion lattice^^ recently discussed 
in conjunction with several ferromagnetic compounds 
MnSi, FeCoSi, and FeGe under magnetic fields jloUii). 
In these compounds, the skyrmion lattice is stabilized by 
the anti-symmetric Dzyaloshinskii-Moriya (DM) interac- 
tion. In contrast, the skyrmion lattice in the triple-(7 
state of the present model is realized via the frustrated 
symmetric exchange interaction. Note that skyrmions of 
our model can take both signs, e.g., skyrmions and anti- 
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FIG. 1. (color online). Phase diagram of the J1-J3 model 
with J1/J3 = —1/3 in the temperature-field plane, obtained 
by a Monte Carlo simulation. Transition temperatures be- 
tween the paramagnetic and the ordered phases are deter- 
mined from the specific-heat-peak position, while those be- 
tween the triple-^ and the single-^ (or the double-^) phases 
are determined by the mixed-phase method ^27.] . 



skyrmions, because our Hamiltonian ([T]) keeps the mirror 
symmetry in the xy spin component. Such Z2 symmetry 
of the Hamiltonian, which is absent in the DM system, is 
spontaneously broken in the triple-g' state. Note that the 
direction of in Eq. (3) is independent to the wavevec- 
tor gr*. In case of the DM system, on the contrary, 
is fixed to be perpendicular to gr* to minimize the DM 
interaction, leading to a helix. 

In the triple-g' state, the state keeps the C3 symmetry 
of the lattice, while it is broken in the single-g' and the 
double-g states. More precisely, Cs symmetry is broken 
only in the xy component in the single-g' state, whereas it 
is broken both in the xy and z components in the double- 
q state. The change in Oi of Eq.(3) induces a translation 
of the skyrmion lattice. 

By comparing the free-energy of these ordered states, 
we construct a mean-field phase diagram in the T-H 
plane. We find that the single-g state always has a lower 
free energy than those of the double-g and the triple- 
q states. Thus, within the mean-field approximation, 
multiple-g states are not stabilized. Meanwhile, we find 
that the free-energy difference between the single-g' and 
the triple-g states becomes small at moderate fields, sim- 
ilarly to the case of a mean-field analysis of the DM ferro- 
magnet [19]. Hence, the fiuctuation effect not taken into 
account in the mean-field approximation might eventu- 
ally stabilize the multiple-g' states. 

In order to further clarify the situation, we perform a 
Monte Carlo simulation based on the standard heat-bath 
method combined with the over-relaxation method. Our 
unit MC step consists of one heat-bath sweep and ten 
over-relaxation sweeps. The lattice is a L x L triangular 
lattice with 36 < L < 288 with periodic boundary con- 
ditions. Typically, a single run contains 2 ~ 4 x 10^ MC 
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FIG. 2. (color online). The intensity plot of the spin struc- 
ture factor for J1/J3 = —1/3, in the single-g' phase (a, b), 
the double-^ phase (c, d), and the triple-g^ phase (e, f). 
The lattice size is L = 72. The upper and lower figures 
represent the spin structure factors of the xy component 
S^{q) = iE^=.,,(IE^^^,Me-^''■'f ) and of the z compo- 
nent S\\[q) = j^{\T.^S^,ze-''^■''^\^). The gray (color) scale 



represents 



steps per spin at each temperature, while averages are 
made over 3^5 independent runs. 

The resulting T-H phase diagram of the Ji- J3 model is 
shown in Fig. [T] for J1/J3 = —1/3. We find that, in addi- 
tion to single-g state, the double-g and the triple-g states 
are stabilized under magnetic fields due to fluctuations. 
In the low temperature limit, the single-g state is always 
stable consistently with the the fixed spin-length condi- 
tion. At larger fields, the single-g phase cuts in between 
the paramagnetic and the double-g' phases. A similar 
phase diagram is obtained also for the J1-J2 model (see 
the online supplement [26]). 

The computed spin structure factors are shown in 
Fig. [2] for each case of the single-g', the double-g and the 
triple-g' phases. As can be seen from the figure, all quali- 
tatively features of the mean- field analysis are met here. 
It should be noticed that sharp spots observed at gr = q* 
are not true Bragg peaks because we are dealing with 
the 2D Heisenberg model. If one recalls that the ordered 
states possess a continuous degeneracy due to the U{1) 
symmetry associated with spin rotations (around z) and 
translations, the observed sharp spots should be quasi- 
Bragg spots associated with power-law spin correlations. 

We show in Fig. [S^a) and (c) typical real-space spin 
configurations obtained in our MC simulation for the 
triple-g state. One can see that the z component, shown 
by the gray (color) scale, often takes an opposite direction 
to the field direction forming a triangular superlattice ex- 
pressed by balck color. The xy component represented 
by arrows forms a vortex (a) or an anti- vortex (c) pattern 
around the black spots. Such a configuration is indeed 
the skyrmion (or anti-skyrmion) lattice predicted by the 
mean- field calculation as a metastable structure (3). In 
Fig. [31(b) and (d), we show the corresponding skyrmion- 
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FIG. 3. (color online). Typical real-space spin configurations 
(left figures) and the intensity plots of the skyrmion density 
(right figures), in the triple-^ phase (a-d) and in the Z phase 
(e, f). Upper (middle) figures represent the skyrmion (anti- 
skyrmion) lattice. The xy component of spins is represented 
by the arrow, while the z component is given by the gray 
(color) scale. Short-time average over 10 MCS is made to 
reduce the thermal noise. The lattice size is L = 72. 



density pattern. The local skyrmion density is defined 
here as the directed area of the sphere surface spanned 
by three spins on every elementary triangle on the lat- 
tice [28]. As can be seen from the figure, in the triple-g' 
state the skyrmion (or anti-skyrmion) forms a triangular- 
lattice with a lattice spacing ^=r-T. The total sum of the 



skyrmion density over the entire system is non-zero, be- 
ing negative for the skyrmion lattice and positive for the 
anti-skyrmion lattice. In the single-g and the double-g' 
phases, this sum turns out to vanish. 

The T-H phase diagram of Fig. [1] contains a new phase 
labeled Z, right to the triple-g' phase, which is not pre- 
dicted in the mean-field analysis even as a metastable 
state. Its existence is suggested from, e.g.^ the specific- 
heat shown in Fig. HJ^a), where clear double peaks are 
observed at 0.39 IJ3I and 0.42 |J3|. The total scalar 

chirality, x = \l {{mT^iXif) with Xi = -^ii • ('^^2 x ^i^) 
where ~ is are three sites on an elementary triangle 
i (both upward and downward), can be regarded as an 
order parameter of the Z2 mirror symmetry. As shown 
in Fig. m^b), x grows at a lower transition temperature, 
indicating that the Z2 mirror symmetry is preserved in 
the Z phase. 

A typical spin configuration in the Z phase is shown 
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FIG. 4. (color online). The temperature dependence of the 
specific heat per spin (a) and of the total scalar chirality per 
plaquette (b) for H/\J^\ — 2.0. The spin structure factor of 
the xy component (c) and of the z component (d) in the Z 
phase. J1/J3 = —1/3 in all cases. 

in Fig. [3l^e). The z component forms a triangular super- 
lattice similar to the triple-(7 phase, whereas the xy com- 
ponent remains disordered (paramagnetic). The corre- 
sponding spin structure factors shown in Fig.|4l^c) and (d) 
also indicate that the xy component exhibits very broad 
spots in contrast to the sharp spots of the z component. 
Thus, in the Z phase only the z component retains a 
quasi-long-range order similar to the triple-g' phase, while 
the xy component remains disordered with a finite corre- 
lation length of about 20 lattice spacings. Skyrmion den- 
sity pattern of the Z phase is shown in Fig. [3jf), where 
one sees that skyrmions and anti-skyrmions are mixed in 
the form of domains. The domains size is comparable to 
the transverse {xy) spin correlation length deduced from 
the width of the peak of S^{q) (Fig. [^c)). Thus, the Z 
phase is a domain state consisting of both skyrmion and 
anti-skyrmion lattices, where the total skyrmion number 
remains zero. At the para-Z transition, the U{1) x U{1) 
symmetry associated with the two independent Oi vari- 
ables of Eq.(3) is broken in an algebraic manner. This 
observation suggests that the para-Z transition is actu- 
ally of the Kosterlitz-Thouless type. 

We now wish to compare the skyrmion-lattice state 
of the present model with that of the DM system as 
observed in MnSi, FeCoSi or FeGe There are 

two important differences between the skyrmion lattice 
of the two systems: (i) First, the skyrmion lattice of the 
present model is much denser than its DM counterpart, 
with smaller lattice constant of ~ 2>/3a (see Fig. [3]). 
In the DM-driven system, the lattice constant is typi- 
cally an order of magnitude larger, because the DM in- 
teraction is usually considerably weaker than the dom- 
inant exchange interaction. Since the anomalous Hall 
conductivity due to skyrmion texture is proportional to 
the skyrmion density [29], larger Hall conductivity is ex- 
pected if the skyrmion (or anti-skyrmion) lattice in the 



triple-g state could be stabilized in appropriate metallic 
materials by the present mechanism, (ii) Second, the Z2 
mirror symmetry is absent in the DM system, while it is 
kept as a Hamiltonian symmetry in the present model, 
being spontaneously broken in the triple-g' state. Its con- 
sequence is that both skyrmion and anti-skyrmion lat- 
tices, inter-connected via the Z2 symmetry, are possible 
in the present model, together with a new Z phase where 
skyrmion and anti-skyrmion lattices form a domain state. 

Triangular-lattice compounds NiGa2S4 and NiBr2 
might be candidates of the multiple-g states. Although 
the observed zero-field properties of NiGa2S4 seem not 
consistent with a simple classical Ji- J3 model in view of 
the absence of an expected first-order transition 
13], its local spin structure is certainly an incommensu- 
rate spiral. Then, the multiple-g' structures might pos- 
sibly be formed under magnetic fields, the required field 
roughly estimated to be 30 50 T. In case of NiBr2, 
estimated value IJ3I/J1 — 0.262 [l^ is much less than 
that of the present study IJ3I/J1 = 3. Although larger 
IJ3I/J1 seems favorable to the formation of the triple- 
q state, further theoretical and experimental studies are 
desirable to clarify the general dependence on | J3 1 / Ji . 

We finally note that within a mean-field approxima- 
tion the multiple-g' states can be obtained as metastable 
states only by assuming the three-fold degeneracy of the 
ordered state. This suggests that the multiple-g' states 
and the skyrmion lattice could be realized not only in 
the triangular lattice, but also in other lattices with a 
trigonal symmetry, e.^., the honeycomb and the kagome 
lattices. 

In summary, we studied the ordering of triangular- 
lattice Heisenberg magnets with an incommensurate spi- 
ral structure. On the basis of a mean-field analysis and 
a Monte Carlo simulation, we found several multiple-g 
phases under magnetic fields in addition to the stan- 
dard single-g phase. The spin structure in the triple-g 
phase is the skyrmion (or anti-skyrmion) lattice. In con- 
trast to the DM-induced skyrmion lattice, the present 
model keeps the Z2 mirror symmetry, which enables 
both skyrmions/anti-skyrmions and gives rise to a new 
Z phase, a domain state consisting of skyrmion and anti- 
skyrmion lattices. 
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FIG. 5. Phase diagram of the Ji — J2 model with J1/J2 = —1/4 in the temperature- field plane, obtained by a Monte 
Carlo simulation. The Hamiltonian of the model is given by Eq. (1) in the main text. Transition temperatures between the 
paramagnetic and the ordered phases are determined from the specific-heat-peak position, while those between the triple-^ and 
the single-g' (or the double-g) phases are determined by the mixed-phase method. 



